Supercritical Coulomb Impurities in Gapped Graphene 
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We study the problem of Coulomb field-induced charging of the ground state in a system of 2D 
massive Dirac particles - gapped graphene. As in its 3D QED counterpart, the critical Coulomb 
coupling is renormalized to higher values, compared to the massless case. We find that in gapped 
graphene a novel supercritical regime is possible, where the screening charge is comparable to the im- 
purity charge, thus leading to suppression of the Coulomb field at nanometer scales. We corroborate 
this with numerical solution of the tight-binding problem on the honeycomb lattice. 
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I. INTRODUCTION 

The successful experimental isolation of a single plane 
of sp 2 carbon atoms (graphene) has stirred our deepest 
assumptions about condensed matter systems 1,2,3 . Not 
only is graphene unique in the realm of solid state, but 
it has also a vast potential as a test ground for many of 
the remarkable predictions of quantum electrodynamics 
(QED) 4,5 . One remarkable characteristic of graphene is 
that its electronic properties can be tailored in many dif- 
ferent ways, either by applying transverse electric and 
magnetic fields, changing its geometry, or modifying the 
substrate where graphene is deposited or grown 3 . 

In fact, while graphene deposited in SiC>2 is well 
described by the two-dimensional (2D) massless Dirac 
equation 1 , graphene grown on SiC can be described in 
terms of massive 2D Dirac electrons 6 . Substrate induced 
potentials can break symmetries of the honeycomb lattice 
and generate gaps in the electronic spectrum. Hence, by 
suitable choice of substrates one can tune the "rest mass" 
of the "relativistic" particles and explore new phenom- 
ena beyond the massless case. Furthermore, the elec- 
tronic properties of devices, such as carrier mobility, de- 
pend strongly on how the electronic degrees of freedom 
(massive or massless) interact with impurities. Of par- 
ticular importance are charged impurities that naturally 
appear either on the substrate, on top of graphene, or 
between the substrate and graphene. Charged impuri- 
ties play an important role in the transport characteris- 
tics of graphene deposited in Si02 7,8,9 and should play 
an important role in epitaxial graphene as well 10 . 

In this paper we explore the non-trivial restructuring 
and charging of the vacuum that occurs in the presence 
of a supercritical Coulomb center 11,12 if the system is 
described by a "massive" spectrum. The problem of an 
unscreened Coulomb impurity in undoped graphene has 
recently received considerable attention 13,14,15,16,17,18 . 
This is justified since, on the one hand, the vanish- 
ing density of states (DOS) at the Fermi energy of un- 
doped graphene and the absence of back-scattering sup- 
press screening 19 . In addition, the Coulomb problem 
in graphene is the condensed matter analogue of su- 
percritical nuclei (Za > 1) in QED, whose rich and 
fundamental phenomena remain elusive to experimental 



testing 12,13,15 . We show here that this analogy achieves 
its fullest in gapped graphene, where one can resolve the 
incremental charging of the vacuum, with strong impli- 
cations for the screening of the Coulomb center. 

The rest of the paper is organized as follows. In Sec. II 
we introduce our model, and in Sec. Ill its properties in 
the massless limit are reviewed. In Sec. IV the behav- 
ior of discrete energy levels in the massive case is dis- 
cussed. We examine in detail the structure of the critical 
wavefunction and the corresponding rcnormalization of 
the critical Coulomb coupling in Sec. V. In Sec. VI we 
study the behavior of the vacuum charge across the criti- 
cal point. Sec. VII contains the corresponding results for 
a finite-size system (where the energy gap is due to the 
finite size only). Sec. VIII contains a discussion of our 
results and conclusions. 



II. THE MODEL 

To address the Coulomb problem in graphene we resort 
to the single- valley effective Dirac description of the elec- 
tron dynamics. We start from the tight-binding Hamilto- 
nian describing electrons in a honeycomb lattice, and in 
the presence of a Coulomb center of strength Z. In addi- 
tion, we assign a different local energy to each sublattice. 
The resulting Hamiltonian is 



H = t 



h.c. 
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In the above, the sums run over unit cells, and the oper- 
ators di(bi) pertain to the A(B) sublattice. This system 
exhibits a band gap of A and, if undoped, has an insulat- 
ing ground state with the Fermi level in the gap. Just as 
in its gapless counterpart, low energy excitations can be 
addressed within an effective mass approximation, con- 
sisting of a k.p expansion around the K and K' points in 
the Brillouin zone. The resulting effective Hamiltonian 
leads to the Dirac equation in 2D under a Coulomb field: 
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where the wavefunction ^(r) has a spinor structure that 
carries the amplitude of the wavefunction on each sub- 
lattice, 



*(r) 



ip B (r) 



(3) 



vf = 3ta/(2h) is the Fermi velocity in graphene and a 
is the C-C distance. It is convenient to introduce g as 
the dimensionless coupling to the external Coulomb field: 
g = Za, and a = e 2 /(hv-p ) is the "fine structure" con- 
stant of graphene. The gap in the spectrum is induced 
by the term proportional to a z , akin to a relativistic rest 
mass. The mass m is related to the local energy mismatch 
in the tight-binding formulation through A = 2mvp . To 
avoid too cumbersome a notation, throughout the paper 
we shall take a system of units in which h — v-p = 1. 
Without loss of generality, the impurity will be consid- 
ered attractive (g > 0) and the C-C distance (a ~ 1.42 A) 
will be used as the length unit. 

In parallel with the exact analytical solution of Eq. (2), 
we solve the tight-binding problem of (1) exactly, via di- 
rect numerical diagonalization on the lattice. This exact 
solution provides an important control of the validity of 
the Dirac approach, and the results of the two approaches 
will be frequently compared. 



III. COULOMB IMPURITIES IN MASSLESS 
GRAPHENE 

Before delving into the details of the supercriti- 
cal regime in the massive case, a brief overview of 
the main physics seen in the massless problem is 
pertinent 13,14,15 . When m = Eq. (2) separates 
in cylindrical coordinates 20 , and can be subsequently 
mapped into the radial Schrodingcr equation of the 
usual Coulomb problem in 3D 15 . One decisive peculiar- 
ity of this mapping is that, although the radial equa- 
tion is formally the same as the radial equation in the 
Schrodinger Coulomb problem, the angular momentum 
quantum number appears replaced by an "effective an- 
gular momentum" equal to I = yj j 2 — g 2 . Here, j is 
the quantum number associated with the total angular 
momentum operator 



1 



(4) 



and takes the values j = ±1/2, ±3/2,. . . . The radial 
solutions are thus expressed in terms of the Coulomb 
functions 15,21 gsign(e), \e\r) and G;(— gsign(e), |e|r). 
Given that, as usual, I determines the asymptotic power 
law behavior of the wavefunctions: 



F,(r~0)~r 



Gi(r~0)- r 



it is evident that g = g c = 1/2 is a singular point for 
the lowest total angular momentum channel (j = ±1/2). 



Below g c the space of solutions is constrained by the re- 
quirement of regularity at the origin, as usual, and this 
selects Fi(r) as the regular solution. But above g c , I 
becomes imaginary and any linear combination of the 
solutions (Fi,Gi) becomes square integrable at the ori- 
gin. At the same time, one sees from the above asymp- 
totics that the solutions oscillate endlessly when r — > 
as cx exp[i log(r)]. This is a signature of the "fall to the 
center" characteristic of highly singular potentials 22 . In 
fact, the supercritical regime can be intuitively under- 
stood from a classical perspective: consideration of the 
classical equations of motion for the Dirac Hamiltonian 
shows that, for a given angular momentum, there will be 
a critical coupling above which the classical orbits spiral 
and fall onto the potential source 13 . The potential has 
become too singular and there are no closed nor scatter- 
ing orbits. 

It is known that, in this case, the quantum mechanical 
problem becomes uniquely defined only after an addi- 
tional boundary condition is introduced, reflecting the 
physical cutoff of the potential at short distances 11,23 . 
For graphene the natural cutoff is the lattice spacing, 
a. This regularization of the potential at short distances 
permits the exact solution to be extended to the super- 
critical regime (g > g c ) which, among other effects, is 
characterized by the presence of marked resonances in the 
spectral density of the hole channel 15 . Such spectral fea- 
tures are analogous to the positron resonances expected 
for a supercritical nucleus in QED 12 . But a fundamen- 
tal difference exists between the two cases: whereas in 
QED, due to the finite mass, the emergence of positron 
resonances is an incremental process (as a function of g), 
in massless graphene an infinite number of them instantly 
appears at g = g c . Scmiclassical considerations illustrate 
how these resonances emerge from an infinite number of 
quasi-bound states embedded in the lower continuum 13 . 
Adding to these spectral peculiarities, the induced elec- 
tronic charge behaves rather differently on the two sides 
of the critical point, being localized close to the impurity 
for g < g c , and otherwise decaying algebraically 13,15 as 
1/r 2 . 

An important question is how will these non- 
perturbative effects contribute to screen the impurity po- 
tential. Adding to the ever present virtual vacuum po- 
larization arising from virtual excitations 18 , in the super- 
critical regime one expects those quasi-localized states to 
partially screen the Coulomb center. It was argued, on 
the basis of a self-consistent treatment, that the impurity 
charge at large distances is screened down to its critical 
value Z c = l/(2a) for any g > g c 13 - We will show that 
in the massive case of interest here the situation is con- 
ceptually different, and a novel regime can be reached 
where the polarization charge is comparable to the bare 
charge Z. This leads to a strong tendency to neutraliza- 
tion of the impurity potential at a length scale set by the 
Compton wavelength in graphene. We examine in detail 
the conditions for this to occur. 
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IV. MASSIVE DIRAC FERMIONS IN 
GRAPHENE 

The emergence of resonant solutions described above 
is best appreciated when the system acquires a mass and 
the electron dispersion becomes e 2 = k 2 + m 2 . This mass 
can arise physically from sublattice symmetry breaking 6 , 
spin-orbit coupling or from reduced dimensionality as in 
a finite-size mesoscopic sample. The one-particle solution 
for the massive case is straightforward and is discussed, 
for instance, in Refs. 14,24. The fundamental difference 
that a finite mass brings to the Coulomb problem is the 
immediate appearance of bound solutions corresponding 
to the hydrogen-like fine structure spectrum in 2D: 



1 sign (g) 



g 2 +[n+ y/j 2 - g 2 



(5) 



In the above n is the principal quantum number and j 
the total angular momentum quantum number defined 
above. These states reside in the gap — m < e < m, 
and are described by wavefunctions that decay expo- 
nentially for distances larger than the "Bohr" radius, 
a o = ^c/{Zoi), where Ac = h/(mv-p) is the "Compton" 
wavelength in graphenc. For an attractive potential the 
lowest bound state (n = 0,j = 1/2) has an energy given 
by 



= mVl-(2ff) 2 , 



(0) 



and hence reaches zero singularly at g = g c , becoming 
imaginary beyond this coupling strength. Similarly to 
what we have discussed regarding the massless case, this 
imaginary energy reflects the fact that a description of a 
point nucleus is impossible for g > g c . A physical regu- 
larization procedure, where one cuts the potential off at 
small r, relaxes this constraint and the bound levels are 
then allowed to sink further 12 ' 25 through negative ener- 
gies until the point e = —m is reached. The diving of 
the bound solutions with increasing coupling is schemat- 
ically depicted in Fig. 1(a). The vicinity of this point, 
where the lowest bound state is about to merge (or has 
just merged) with the continuum is the focus of all our 
subsequent discussion. 



V. CRITICAL WAVEFUNCTION AND 
CRITICAL COUPLING RENORMALIZATION 

Given that in a regularized Coulomb potential the 
bound levels are allowed to dive into negative energies, 
the relevance is transferred from the particular coupling 
g = g c = 1/2 to the value g c at which eo = —m. Nec- 
essarily g c > g c , and the critical coupling is thus renor- 
malized to higher values. This is one of the immedi- 
ate consequences of the presence of a gap. The merging 
of the bound solution with the continuum is nonethe- 
less peculiar, and hence is worth carcfull consideration. 
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FIG. 1: (a) Schematic depiction of the change of bound levels 
(5) with increasing coupling. The dashed line represents the 
evolution of eo for a point nucleus, and the solid one is for 
a regularized potential. Adapted from Ref. 26. (b) Effective 
semiclassical momentum of Eq. (8) when eo ~ — m and g > g c 
(solid) or g < g c (dashed). 



We begin with a semiclassical argument. The first im- 
portant thing to notice is that the merging state does 
not become completely delocalized, as one would expect 
in typical Schrbdinger problems, but remains localized 
as it dives into the lower continuum. This is somewhat 
counter-intuitive and is best appreciated by considering 
the problem scmiclassically within WKB. An alternative 
to the full WKB approximation is to consider the "clas- 
sical relativistic" momentum 



p{r) 2 = [e-U{r)] 2 - m 2 



(7) 



where U (r) is the Coulomb potential. With this defini- 
tion, the radial momentum p r is written as 



p r {r) 2 = k 2 + 2ge/r + (g 2 - i 2 )/r 2 



(8) 



with 

,2 



the classical effective angular momentum and k 2 = 
e" — m 2 . The solid line in Fig. 1(b) shows the profile of 
p r {r) 2 slightly after the merging has taken place. There is 
a classically forbidden region [p(r) 2 < 0] between r_ and 
r + , and the long range Coulomb tail implies that r + 3> 
r_. In fact, precisely at eo = — m we have r+ = oo, and 
the barrier is infinite. Close to the critical point we can 
write eo — — m — (3(g — g c ) and examine the penetrability 
of this barrier within WKB. It reads 



e 25cl = exp 



2-Kmg c 



y/2pm*fg - g c 



(9) 



and becomes exponentially small as g — > g£. Such wide 
barrier (absent when eo — +m) ensures that the state 
remains appreciably localized, even after merging into 
the lower continuum. 

We now turn to the exact quantum mechanical solution 
of Eq. (2) for the particular regularization of the Coulomb 
potential described by 25 



V(r) 



-g/r, r>a 
-g/a, r<a 



(10) 
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This regularization represents the physical situation in 
which there is a Coulomb impurity at the center of an 
hexagon in the honeycomb lattice. In this case the lat- 
tice parameter, a, is the closest distance that an electron 
hopping among carbon sites can be from the potential 
source. It is also a good approximation for an impurity 
placed slightly above the graphene plane. Cylindrical 
symmetry is preserved by this regularization and Eq. (2) 
naturally separates in cylindrical coordinates. We define 
the radial spinor components as (j = ±1/2, ±3/2, . . . ) 



*.fr«rt- 1 ( e - iU - 1/2)VA ( r ) 



(11) 



after which the radial equations for (2) become (r > a) 

0. (12) 



(e+f-m) + 
(dr-i) (e+f+m). 



A(r) 
B{r) 



We will be interested in the states at the threshold, with 
energy e = — m. In that case, the equation for the A{r) 
component reads (r > a) 



r 2 (^A(r))" 



2 -2 , 
9 -J + 



1 



2gmr 



y/rA{r) = , (13) 



whose solution is given in terms of the modified Bessel 
functions I v (z) and K v (z). By imposing a vanishing 
boundary condition at infinity one obtains (Af is a nor- 
malization constant), 



A(r) =MK iu ( y /%gmr), 



(14a) 



where v = 2y 'g 2 — p. We note that this solution has 
the same form as the corresponding spinor component 
in the 3D QED problem 11 - 25 . The B component follows 
directly from (12): 



B(r) 
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3 + l -^)K iv y%gmr) + 



y / 2gmrK tu _ i ( a/ 8gmr) 



(14b) 



The solutions (14) pertain to the region r > a, where 
the potential (10) has the Coulomb form. For distances 
smaller than the regularization distance the solution for 
the radial spinor components of (11) is given in terms of 
cylindrical waves: 



A(r) 
B(r) 



gJj-i/2(kr) 
kaJ j+1 / 2 (kr) 



(15) 



where we have introduced ka = \J g 2 — 2mga. The solu- 
tions (14) and (15) are at energy e = — m and need to be 
matched at r 



(16) 



= a: 
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FIG. 2: (color online) Critical coupling, g c , as a function of 
the mass/gap obtained from the solution of Eq. (16). The top 
horizontal axis shows ma in the units natural for the tight- 
binding calculation: ma — > a/Ac = A/(3t). The inset shows 
the probability density associated with the critical wavefunc- 
tion at e = — m [cfr. Eq. (17)]. 



This completes the determination of the critical spinor 
wavefunction, 'l'c(f)- Since the Dirac equation (2) is 
of first order in the gradient operator, one is required 
to match only the spinor components (and not, addi- 
tionally, their derivatives as would be the case for the 
Schrodinger equation), and this is enough to guarantee 
the continuity of the probability current density across 
the region r = a. The matching procedure generates a 
transcendental equation for <?, with an infinite number 
of solutions, on account of the oscillating character of 
the functions Ki U (z). For each angular momentum j, 
the multiple solutions correspond to the values of g for 
which the higher bound states (i.e. those levels having 
the same j-symmctry) reach the continuum. To deter- 
mine g c we concentrate on the s level (j = 1/2), and 
solve Eq. (16). Its smallest solution for g yields g c , the 
renormalized critical coupling. A plot of g c as a function 
of ma is shown in Fig. 2. It can be seen that the depar- 
ture of g c from the value g c = 1/2 is singular at the origin, 
which means that even an arbitrarily small mass causes 
a significant change in the critical coupling. In the figure 
we also point out (dashed lines) the value of g c ~ 0.949, 
that corresponds to a gapped graphene spectrum with 
A = 2mvp = 0.2t (in units of the tight-binding hopping 
parameter). We have chosen this value for illustration 
purposes, to be discussed later. The critical probability 
amplitude, defined in terms of the critical wavefunction 
as p(r) — ^l(r)^> c (r) , is shown in the inset. The crux of 
our argument resides in the fact that ^ c is clearly local- 
ized, as can be inspected from the asymptotic behavior 
of p(r) for r 3> 1/m: 



p(r) = *J(r)¥ c (r) 



-1/2 

r ' exp 



mg c r 



(17) 
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Restoring the units, the "Compton" wavelength emerges 
as the characteristic localization length: mr — > 
mVF r/h = r/Xc- 

In order to ascertain the validity of these results in the 
context of the original tight-binding problem, we explic- 
itly compare the above with the exact numerical solution 
in the honeycomb lattice. Of particular interest are the 
renormalization of the critical coupling and the character 
of the diving states. In Fig. 3 we present the exact tight- 
binding spectrum for a Coulomb center in the honeycomb 
lattice with a sublattice energy mismatch (energy gap) 
of A = 0.2t. Inspection of the main panel reveals sev- 
eral features, among which: (i) the lowest positive level 
(which corresponds to eo) immediately detaches from the 
continuum and sinks rapidly, followed by other states at 
higher g; (ii) the level eo reaches — m at g = 0.83 and 
touches the lower continuum at g = 0.87, in a very good 
agreement with the result g c ~ 0.9 predicted in Fig. 2 
from the solution of the Dirac equation 27 ; (iii) as g in- 
creases past <? c this level sinks further, as is evident from 
the level avoidances highlighted by the dashed/shaded re- 
gion; (iv) the critical wavefunction, shown in the inset, is 
highly concentrated around the origin, as expected from 
the foregoing, namely Eq. (17). The lowest bound level 
in Fig. 3 is doubly degenerate. This is expected since 
although the lowest bound solution of the Dirac equa- 
tion (6) is non-degenerate, there is a valley degeneracy 
to account for in the Dirac description. 

While the continuum solutions are sensitive to the fi- 
nite size of the numerical system 27 , one does not expect 
the lowest bound solution to be markedly sensitive for the 
sizes used in this calculation. Hence we can extrapolate 
the trajectory of the lowest bound state in Fig. 3 to the 
thermodynamic limit. In that case, the critical coupling 
in the lattice would be at ~ 0.83, which is smaller than 
the value g c = 0.949, and indicates that the effective cut- 
off distance for the regularization (10) is slightly smaller 
than a. 



VI. VACUUM POLARIZATION AND VACUUM 
CHARGING 

It is clear, either from Eq. (17) or from the exact re- 
sults in the lattice [Fig. 3 (inset)], that the merging state 
has a localized character. The characteristic length scale 
Ac is related to the gap through Ac /a = 3t/A. For the 
gap A = 0.26 eV, reported in reference 6 this corresponds 
to Ac — 30a. However, the size of the merging bound 
state, measured as the average radius and calculated by 
using the critical wave-functions, is smaller: (r) = 13.9a. 
Beyond this scale, the impurity potential is screened by 
essentially one charge unit, times the product of the spin 
and valley degeneracies (N = 4) . Of course this screening 
is meaningfull only in a very diluted impurity configura- 
tion (rij <C 10 12 cm~ 2 for the quoted experimental gap). 

We can appreciate this more clearly by studying the 
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FIG. 3: (color online) Low energy close-up of the exact energy 
spectrum (E n ) as a function of the coupling g for the tight- 
binding problem (1). The inset contains the exact numerical 
wavefunction ^E , c ('"), with the dashed line marking Ac for this 
gap. The sublattice gap is A = 0.2t, and a lattice of 124 2 
sites with a central impurity has been used. 

induced charge, defined as 
Sp(r)= E X ] E {r)XE{r)- £ X°e (r) X ° E (r) , (18) 

E<E F E<-m 

where XE(f) are the continuum wavefunctions in the 
presence of the potential, and X E ( r ) stand for the contin- 
uum wavefunctions in the absence of potential. We work 
with a constant number of electrons as the potential is 
turned on, and hence Ep ^ —m in general 28 . We need 
also to consider the intermediate situation in which there 
is a bound state just about to merge with the lower band. 
The continuum wavefunctions in this case are labeled as 
XbM- F° r the sake of the argument assume that there is 
only one bound state that we denote by f c (r) and, fur- 
thermore, let us accept that we can project everything 
onto the states E < m. In other words, we accept that 
the states E < m approximately form a complete set in 
each circumstance: 

E X% f (r)x° E (r')^S(r-r') (19a) 

E<-m 

E XEHr) X E(r')^5(r-r') (19b) 

E<-m 

E XE\r)x E {r') + *J{r)* c {r')~8{r-r'). (19c) 

E<-m 

From these relations follows that, above g c when one 
state has dived onto the continuum, the induced charge 
at constant density defined in (18) can be approximated 
as 

Sp(r) ~ |* c (r)| 2 + Sp pol (r) - | X - m (r)| 2 , (20) 
where Sp po i(r) represents the polarization charge caused 
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by the deformation of the continuum wavefunctions only: 
S Ppol (r)= E \x c E (r)\ 2 -\x%(r)\ 2 , (21) 

E<-m 

with 

J Sp po i(r)dr = 0, 

and X- m ( r ) represents the topmost plane wave, that ap- 
pears because we keep the number of electrons constant 
as g is varied. Expression (20) is valid for one bound 
level merging into the continuum but can be generalized 
for the additional subsequent divings. 

The result (20) shows that, although after the merg- 
ing there is (strictly) no localized state, the total induced 
charge (20) remains essentially determined by the profile 
of the critical state, ^> c (r), plus the background polariza- 
tion charge. The last term in eq. (20) is a phase-shifted 
wave, and thus the smallest and neglcctable contribution 
to Sp(r). The constant number of electrons, and the fact 
that all states are normalized, evidently implies that the 
integral of Sp(r) over the entire volume vanishes. This 
allows for the definition of Q{R), the total charge inside 
a radius r = R: 

Q(R) = N [ Sp(r)dr. (22) 

J\r\<R 

Here N represents the degeneracies of the problem: N = 
4 (spin x valley) for the solution of the single- valley Dirac 
equation, and N = 2 (spin) for the numerical solution 
in the honeycomb lattice. This function Q(R) will raise 
quickly at small R, rapidly attaining a maximum on ac- 
count of the localized profile of ^ c (r) 29 . Q(R) then re- 
turns slowly to zero as R is increased further, because 
of the normalization of the states and the fact that we 
work at a constant number of electrons. The maximum 
in Q(R) (which we designate Q m ax) is thus a measure of 
the amount of charge pulled to the vicinity of the impu- 
rity, and is our parameter of interest. For g < g c the 
only contribution to Q(R) comes from the polarization 
charge, Sp po i(r); but whenever a new level dives, a con- 
tribution of the type |\E' c (r)| 2 adds to the total induced 
charge and one should observe discontinuous jumps in 
Qmax with increasing g. 

Once again this expectation is confronted with exact 
results on the lattice. We compute (18) numerically using 
the exact wavefunctions of the tight-binding model. In 
Fig. 4(a) we plot the resulting Qmax- Up to g ~ 0.9 the 
only contribution comes from Sp po i and the variation of 
Qmax is smooth. At the point where in Fig. 3 the first 
level merges with the continuum the first discontinuity 
occurs in Q m ax, which jumps by roughly 4 units (N x 1 = 
4), just as expected. Also as expected, this charge is 
concentrated within a region of radius ~ Ac- Therefore, 
the approximations used to obtain expression (20) are 
legitimate, and the induced charge beyond g c is mostly 
determined by the profile of the critical states. 
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(b) 

FIG. 4: (color online) (a) Plot of Q ma x (the amount of charge 
pulled to the vicinity of the impurity, see text) versus g from 
exact numerical diagonalization of the tight-binding problem 
[Fig. 3]. The steps signal the incremental diving of bound 
levels into the lower continuum. The distances at which 
Q(R) — Qmax are indicated for the first two steps, (b) The 
same finite system without an explicit gap (A = 0). In 
both cases the system has 124 x 124 carbon sites and lin- 
ear dimensions of 107a x 186a. The insets amplify the low 
g region, which is compared with the RPA result for m—0: 
Q(R) = f g J ]<R 5(r)dr = § g (dashed line). 



This is the precise analogue of the QED prediction 
for the charging of the vacuum. In a superheavy nu- 
cleus with high enough Z the hydrogen-like bound lev- 
els can merge with the positron continuum. When that 
happens, the Schwinger mechanism of electron-positron 
creation becomes spontaneous since the bound level has 
become degenerate with the positron states, and conse- 
quently there is no energy cost involved in creating an 
electron-positron pair. The electron occupies the bound 
level close to the nucleus, while the positron escapes to 
infinity 12 . Since the positron continuum constitutes the 
vacuum of QED, the spontaneous creation of an electron 
in a level below — m leads to a restructuring of the vac- 
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uum that acquires a charge Q = — 2e 30 . The curve in 
Fig. 4(a) mimics entirely the QED prediction found, e.g., 
in Fig. 1.5 of Rcf. 12. The striking difference rests in 
the magnitude of the effect. In QED, due to the small- 
ness of the real fine structure constant (q;qed — 1/137) 
the nuclear charge required to reach the critical regime is 
very large (Z c ~ 170), and the jumps in the polarization 
charge lead to an effective atomic number Z cS = Z c — Q, 
with Q <C Z c : a small correction. 

In graphene, for estimate purposes, we assume a SiC 
substrate (with dielectric constant e ~ 10), as used in 
Ref. 6, leading to a ~ 0.4. For the gap A = 0.26 eV 
seen in this experiment we obtain g c = (Za) c = 0.84 
and, consequently, the critical valence is Z c ~ 2.1. This 
means that, according to our discussion, impurities with 
valence Z > 2 are expected to be completely screened 
beyond the scale (r) = 13.9a « 2nm, since at this dis- 
tance the amount of accumulated charge in the vacuum 
is Q = 4, and thus Z cS (r > (r)) = Z-Q « 0. In fact, for 
Z < 4 over-screening takes place, and we expect many- 
body interactions to be important, effectively removing 
the over-screening tendency; however such calculations 
are very complex and beyond the scope of this work. In 
addition, we note that substrates with smaller dielectric 
constants lead to higher values of Za and are, in princi- 
ple, capable of reducing the valences Z needed to observe 
the effect to Z = 1,2. 

Electron-electron interactions can also lead to a change 
of the critical coupling. Estimated perturbativcly (a <C 
1) this change is Sg c ~ ag c , which would not affect sig- 
nificantly the physics discussed here (but would clearly 
affect the precise value of Z c for a given a). We recall 
also that the above theory assumes the chemical poten- 
tial to be at n = —m, while the experiments of Ref. 6 in 
fact have /i > +m. Tuning the Fermi level to the gap via 
chemical doping or gating and varying a by changing the 
substrate's dielectric constant would provide a possible 
way of exploring our predictions experimentally 



VII. VACUUM CHARGING IN A FINITE 
SYSTEM 

The discrete eigenstates of a finite graphene sheet, and 
the fact that the density of states vanishes at the Dirac 
point, lead to another realization of a gapped spectrum, 
and the question of whether the effects discussed above 
carry to this case arises naturally. We have diagonalized 
the same lattices as above, this time without any sublat- 
tice mismatch (i.e.: A = 0, or m — 0), corresponding to 
the massless Dirac problem on a finite system. In this 
case we found that the gap induced by the finite size 
of the system is A ~ 0.06t, Ac ~ 50a, and (r) « 21a. 
This leads to a renormalizcd g c ~ 0.76, representing a 
50% increase with respect to g c just from the fact that, 
although massless, the system is constrained to a finite 
size. For brevity we show only the integrated polariza- 
tion charge Q max in Fig. 4(b). The curve of <5 ma x displays 



the characteristic step discontinuity at 0.7 < g < 0.8, in 
agreement with the estimate for g c . Of course, in this 
case the situation is peculiar since the effective Ac is tied 
to the linear size of the system, L, and the gap is now 
A cx Consequently Ac, although still a fraction of 

L, grows with the system size. Nevertheless, most impor- 
tantly, the picture of vacuum charging discussed above 
remains valid, with the vacuum charge residing at a dis- 
tance (r) « 21a, significantly smaller than the linear size 
L ~ 107a. Thus the vacuum charging can be potentially 
observed in mesoscopic samples as well. 



VIII. DISCUSSION AND CONCLUSIONS 

The fact that gapped graphene is described in terms 
of massive Dirac quasiparticles makes it a rather unique 
solid state system. We have seen, in particular, that ef- 
fects of supercritical vacuum charging are expected when 
undoped graphene is exposed to low-valence charged im- 
purities. 

One of the key steps in our calculation is the explicit 
regularization of the potential (10) that allows the div- 
ing of the bound levels into the lower continuum. Since 
an impurity located slightly out of the plane would still 
lead to a regularized Coulomb potential, our results are 
not altered qualitatively because details of the regular- 
ized potential at short distances are not important 11 . At 
the quantitative level, the regularization distance deter- 
mines the renormalization of g c . Since the typical dis- 
tances to the plane expected for chemisorbed impurities 
in graphene are still in the sub- nanometer range 9 , we do 
not expect a prohibitive increase in g c nor, consequently, 
in Z c . Therefore the charged impurity does not need to 
be embedded in the graphene plane, and can be simply 
adsorbed to its surface. 

The theory predicts a strong tendency towards screen- 
ing of the external Coulomb potential on nanometer 
scales, the actual values depending on the details of the 
particular situation (notably, the magnitude of the gap). 
We mention also that, despite our initial assumption of 
attractive impurities (g > 0), this is only for convenience. 
The particle-hole symmetry of the problem ensures that 
the results remain valid for repulsive charges (negative 
ions), in which case the (positron) bound levels merge 
with the upper continuum of states at precisely the same 
critical couplings. 

Whereas in gapless graphene the supercritical regime 
is characterized by an infinite number of resonances in 
the hole (positron) channel 13,15 , the phenomenon of vac- 
uum polarization in gapped graphene is an incremental 
process. Each level dives at successively higher values of 
the coupling g = Za, leading to the the step- wise shape 
of the curves of Q max shown in Fig. 4. This translates 
to a quite dissimilar screening of the supercritical impu- 
rity in the massless and massive cases: for the latter we 
can have complete screening at very short distances in a 
system with no carriers. 
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Finally, the estimates of the critical valence Z c ~ 2 
made in this paper are based on the specific parameters 
of the epitaxial samples used in Ref. 6 (to wit the gap A 
and the dielectric constant of the substrate) . The super- 
critical regime is determined by Za > g c , to which several 
players contribute: Z itself, the dielectric medium (via 
a), the mass/gap and the regularization distance (via g c ). 
Some of these are more manageable to experimental con- 
trol than others. Encouraging experiments are emerging 
that seek control over some of them. In Ref. 9 controlled 
doping with monovalent ions has been achieved, and pre- 
sumably the same could be done with divalent alkaline 
ions. In Ref. 31 the fine structure constant a in exfoli- 
ated graphene could be controlled through changes in the 
dielectric environment. And in Ref. 32 a gap was found 
for graphene flakes on Ni surfaces. 

The charging of the vacuum in the presence of a strong 
Coulomb center is a long standing prediction of QED in 



strong fields that remains unconfirmed through a direct 
experiment. This is an obvious consequence of the diffi- 
culties imposed by the nuclear charge of Z ~ 170 required 
to reach the supercritical regime in QED. Our calcula- 
tions suggest that, under the conditions discussed, the 
analogous effect might be observable in a solid-state con- 
text, with low-valence ions. 



Acknowledgments 

We are grateful to B. Uchoa, O. Sushkov, and R. 
Barankov for stimulating discussions. V.M.P. is sup- 
ported by FCT via SFRH/BPD/27182/2006 and POCI 
2010 via PTDC/FIS/64404/2006, and acknowledges 
Centro de Ffsica do Porto for computational support. 



1 K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, 
M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and 
A. A. Firsov, Nature 438, 197 (2005). 

2 A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. 
Novoselov, and A. K. Geim, arXiv:0709.1163 (2007). 

3 A. K. Geim and K. S. Novoselov, Nat. Mater. 6, 183 (2007). 

4 A. H. Castro Neto, F. Guinea, and N. M. R. Peres, Physics 
World 19, 33 (2006). 

5 M. I. Katsnelson and K. S. Novoselov, Solid State Com- 
mun. 143, 3 (2007). 

6 S. Y. Zhou, G.-H. Gweon, A. V. Fedorov, P. N. First, W. A. 
de Heer, D.-H. Lee, F. Guinea, A. H. Castro Neto, and 
A. Lanzara, Nature Materials 6, 770 (2007). 

7 K. Nomura and A. MacDonald, Phys. Rev. Lett. 98, 
076602 (2006). 

8 S. Adam, E. H. Hwang, V. M. Galitski, and S. D. Sarma, 
Proc. Nat. Acad. Sci. 104, 18392 (2007). 

9 J. H. Chen, C. Jang, M. S. Fuhrer, E. D. Williams, and 
M. Ishigami, Nature Physics 4, 377 (2008). 

10 W. A. de Heer, C. Berger, X. Wu, P. N. First, E. H. Con- 
rad, X. Li, T. Li, M. Sprinkle, J. Hass, M. L. Sadowski, 
et al., Solid State Comm. 143, 92 (2007). 

11 Y. B. Zeldovich and V. S. Popov, Sov. Phys. Usp. 14, 673 
(1972). 

12 W. Greiner, B. Muller, and J. Rafelski, Quantum Electro- 
dynamics of Strong Fields (Springer, 1985), 2nd ed. 

13 A. V. Shytov, M. I. Katsnelson, and L. S. Levitov, Phys. 
Rev. Lett. 99, 236801, 246802 (2007). 

14 D. S. Novikov, Phys. Rev. B 76, 245435 (2007). 

15 V. M. Pereira, J. Nilsson, and A. H. Castro Neto, Phys. 
Rev. Lett. 99, 166802 (2007). 

16 R. R. Biswas, S. Sachdev, and D. T. Son, Phys. Rev. B 
76, 205122 (2007). 

17 M. M. Fogler, D. S. Novikov, and B. I. Shklovskii, Phys. 
Rev. B 76, 233402 (2007). 



18 I. S. Terekhov, A. I. Milstein, V. N. Kotov, and O. P. 
Sushkov, Phys. Rev. Lett. 100, 076803 (2008). 

19 T. Ando, J. Phys. Soc. Jpn. 75, 074716 (2006). 

20 D. P. DiVincenzo and E. J. Mele, Phys. Rev. B 29, 1685 
(1984). 

21 M. Abramowitz and I. A. Stegun, Handbook of Mathemat- 
ical Functions (Dover, New York, 1964). 

22 L. D. Landau and E. M. Lifshitz, Quantum Mechanics: 
Non-Relativistic Theory (Pergamon Press, 1981). 

23 K. M. Case, Phys. Rev. 80, 797 (1960). 

24 V. R. Khalilov and C.-L. Ho, Mod. Phys. Lett. A 13, 615 
(1998). 

25 V. S. Popov, Sov. J. Nucl. Phys. 12, 235 (1971). 

26 B. Muller and J. Rafelski, Phys. Rev. Lett. 34, 349 (1975). 

27 In Fig. 3 the boundaries of the continuum also change with 
g. This is due to the fact that, in such a finite system as the 
one analyzed there, all the energy levels are phase-shifted 
by the potential. This remains true for the levels close to 
the gap and therefore the merging of the bound solution 
with the continuum in a finite-sized lattice appears at a 
larger coupling. This phase shift effect vanishes when the 
size of the system becomes infinite. 

28 The chemical potential (the last occupied level) follows the 
bold black line in Fig. 3. 

29 There is an additional localized contribution that comes 
from Sppd which is reminiscent of the localized screening 
charge oc 8(r) in the massless case 18 . This term yields the 
smoothly varying part in Fig. 4, but is not our focus here. 

30 The factor of 2 comes from the spin degeneracy. 

31 C. Jang, S. Adam, J.-H. Chen, E. D. Williams, S. D. 
Sarma, and M. S. Fuhrer, arXiv:0805.3780 (2008). 

32 A. Gruneis and D. V. Vyalikh, Phys. Rev. B 77, 193401 
(2008). 



